/* Kiefer.c - simulates two Kiefer processes and extracts
   a certain functional of it. */

#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <sys/time.h>
#include <time.h>

#define N 50 /* The number of simulations preformed*/
#define dt 60 /* The grid-resolution 1/dt uniformly */
#define M 30 /* Number of elements in the sum */
#define PI M_PI

int main(int argc, char* argv[])
{

  int times = 0;

  int t1 = 0;
  int t2 = 0;
  int s = 0;
  
  int i1 = 0;
  int i2 = 0;
  int k = 0;

  double result1 = 0;
  double result2 = 0;
  double constant = 0;

  double static diffK[dt][dt][dt]; 
  double static normal1[M][M][M];
  double static normal2[M][M][M];
  double simulations[N];

  double maxK[dt];
  double maxtotal = 0;
  
  FILE * pFile; /* File pointer for saving our results */

  /* Gsl variables: */
  gsl_rng *r;
  double t;
  struct timeval tim;
  
  /* Initialize the random number generator
     with a seed given by the time at the start of the 
     simulation */
  gettimeofday(&tim, NULL);
  t=tim.tv_sec+(tim.tv_usec/1000000.0);
  
  /* We use the MT19937 random number generator
     by Makoto Matsumoto and Takukji Nishimura */
  r=gsl_rng_alloc(gsl_rng_mt19937);
  gsl_rng_set(r, t);

  printf("Computation started. \n");
  
  for(times = 0; times < N; times++) {

    /* First, we generate the array of normal
       variables to use in the expansion process */
    for (i1 = 1; i1 < M; i1++) {
      for (i2 = 1; i2 < M; i2++) {
	for (k = 1; k < M; k++) {
	  normal1[i1][i2][k] = gsl_ran_gaussian(r,1);
	  normal2[i1][i2][k] = gsl_ran_gaussian(r,1);
	  
	}
      }
    }
    /* We now apply the Karhunen-Loeve expansion */

    for (t1 = 1; t1 < dt; t1++) {
      for (t2 = 1; t2 < dt; t2++) {
	for (s = 1; s < dt; s++) {
	  result1 = 0;
	  result2 = 0;
	  for (i1 = 1; i1 < M; i1++) {
	    for (i2 = 1; i2 < M; i2++) {
	      for (k = 1; k < M; k++) {
		constant = (M_SQRT2 * sin(((double) i1)* M_PI* 
		  (((double)t1)/dt))/(((double)i1) * M_PI)) 
		  * (sqrt(2) * sin( ((double)i2)* M_PI* 
		  (((double)t2)/dt))/(((double)i2) * PI))  
		  * sqrt(2) * sin((((double)k) + 0.5)*M_PI
                  *(((double)s)/dt))/((((double)k) + 0.5)*M_PI);
		result1 = result1 + 
		  normal1[i1][i2][k] *constant;
		result2 = result2 + 
		   normal2[i1][i2][k] *constant;
	      }
	    }
	  }

	  /* We discard the specific path realiaztions 
	     of the Kiefer process as it is not what is
	     of interest to us in the current situation.
	     This would not be the case in other settings,
	     and this is the place to save the paths of
	     our simulations */
	  diffK[t1][t2][s] 
	    = fabs(result1 - result2);
	}
      }
    }

    /* We now have two realizations of a Kiefer-process
       and continue to extract a certain functional 
       of them of K.     
       In specific, we want to find 
       ( \sup_{0 \leq s \leq 1}  \sup_{t\in [0,1]^2} 
       | K1(s,t) - K2(s,t)|)^2 
       for this approximation of K1 and K2.
       Thus, we first extract the first supremum
       for each s */
    for (s = 0; s<dt; s++) {

      /* note that diffK[s][0][0] is not initialized,
	 and would in any case be zero from the structure
	 of Kiefer-process.*/
      maxK[s] = diffK[s][1][1];
      
      for(t1 = 1; t1 < dt; t1++) {
	for(t2 = 1; t2 < dt; t2++) {
	  
	  if (diffK[s][t1][t2] > 
	      maxK[s]) {
	    maxK[s] = 
	      diffK[s][t1][t2];
	  }	  
	}
      }
    }
    /*We now extract the supremum over s */ 
    maxtotal = maxK[0];
    for (s = 0; s<dt; s++) {
      if(maxK[s] > maxtotal) {
	maxtotal = maxK[s];
      }
    }
    /* Finaly, we are interested in
       the squared maximum! */    
    simulations[times] = maxtotal*maxtotal;
    
    printf ("\n simulations[%d] = %g\n", times, 
	    simulations[times]);  
  }

  printf("Computation finished.\n");

  /* Saves the result. */
  pFile = fopen(argv[1], "w+");
  for(times = 0; times < N; times++) {
    if (times < N-1) {
      fprintf(pFile, "%f, ", simulations[times]);
    }
    else {
      fprintf(pFile, "%f \n", simulations[times]);
    }
  }
  fclose(pFile);
  
  gsl_rng_free(r);
  
  return 0;
}
